# -*- coding: utf-8 -*-
"""
@author:xuyuntao
@time:2021/1/13:10:47
@email:xuyuntao@189.cn
"""
import numpy
import numpy.fft as nfft
import scipy.signal as scisignal
import csv

dataPath="AnalogModulationSignal"
for fileNum in range(1):
    data=[]
    with open(dataPath+"\\"+"data{}.csv".format(fileNum+1),"r") as f:
        csvRead=csv.reader(f)
        for _ in csvRead:
            data.append(_)
    signal=numpy.array(data,dtype=numpy.float32).reshape(-1)

    print("--------文件：{}--------".format("data{}.csv".format(fileNum+1)))

    # ---------常量---------
    carrier_f = 10e6  # 载波频率，指经过混频后的中频信号，10MHz
    f_sample = 8 * carrier_f  # 采样频率，8倍的载波频率
    T_sample = 1 / f_sample  # 采样周期
    # signalLength = signal.size
    t_a = 1
    t_gamma_max = 500 # 5.8  # gamma_max的特征门限值
    t_sigma_dp = 0.5  # sigma_dp的特征门限值
    t_P = 0.6  # P的特征门限值
    n_series=numpy.linspace(0,signal.size-1,signal.size)
    time_series=n_series*T_sample
    Amp_pass = 1 # Ap：通带最大衰减
    Amp_stop = 30 # As：阻带最小衰减

    # ---------载频测量---------
    # 基于过零点检测估计频率
    signal_augment=numpy.zeros(signal.size+1)
    signal_augment[1:signal.size+1]=signal
    signal_forword = signal_augment.copy()
    signal_augment[0:signal.size]=signal
    signal_augment[-1]=0
    signal_after = signal_augment.copy()
    zero_crossing = (signal_forword*signal_after) < 0
    zero_crossing_num = numpy.sum(zero_crossing[2:signal.size])

    # 计算载频
    last_zero_2 = numpy.where(zero_crossing[0:signal.size])[0][-1]
    last_zero_2_value = signal[last_zero_2]
    last_zero_1 = last_zero_2 - 1
    last_zero_1_value = signal[last_zero_1]

    first_zero_2 = numpy.where(zero_crossing[0:signal.size])[0][0]
    first_zero_2_value = signal[first_zero_2]
    first_zero_1 = first_zero_2 - 1
    first_zero_1_value = signal[first_zero_1]
    z_last = last_zero_2 - last_zero_2_value * (last_zero_1_value - last_zero_2_value) # 最后一个过零点区域斜率
    z_first = first_zero_2 - first_zero_2_value * (first_zero_1_value - first_zero_2_value)
    f_signal = 0.5 * f_sample * (zero_crossing_num - 1) / (z_last - z_first)
    f_Res = f_sample / signal.size # 频率分辨率
    kc = round(f_signal / f_Res) # 中心频率f0在X(k)中的位置

    # ---------计算P - 谱对称性---------
    signal_fft = nfft.fft(signal)
    P_L = numpy.sum(numpy.abs(signal_fft[0:kc])** 2)
    P_U = numpy.sum(numpy.abs(signal_fft[kc+1:2*kc+1])** 2)
    P = (P_L - P_U) / (P_L + P_U)

    # ---------正交解调---------
    w_signal = 2 * numpy.pi * f_signal # 模拟角频率
    signal_inPhase = signal* 2* numpy.cos(w_signal*time_series)
    signal_quadPhase = signal* -2*numpy.sin(w_signal*time_series)

    # ---------低通滤波---------
    f_pass = round(f_signal) # fp：滤波器的通带频率
    f_stop = round(0.9 * f_signal) # fst: 滤波器的截止频率
    f_pass_digital = 2 * f_pass / f_sample
    f_stop_digital = 2 * f_stop / f_sample
    filterSteps,omega_n=scisignal.buttord(f_pass_digital,f_stop_digital,Amp_pass,Amp_stop)
    filter_b,filter_a=scisignal.butter(filterSteps,omega_n,"low",output="ba")
    signal_inPhase_lowPass=scisignal.filtfilt(filter_b,filter_a,signal_inPhase)
    signal_quadPhase_lowPass=scisignal.filtfilt(filter_b,filter_a,signal_quadPhase)

    # ---------计算各个特征值---------
    phi=numpy.arctan2(signal_quadPhase_lowPass,signal_inPhase_lowPass)
    phi_avg = numpy.mean(phi) # 平均值
    phi_NL = phi - phi_avg
    amp = numpy.sqrt(signal_inPhase_lowPass**2 + signal_quadPhase_lowPass**2) # 幅度
    amp_avg = numpy.mean(amp) # 平均值
    amp_n = amp / amp_avg # 用平均值对瞬时幅度进行归一化
    amp_cn = amp_n - 1

    # ---------计算特征值gamma_max：归一化零中心瞬时幅度之谱密度的最大值---------
    gamma_max=numpy.max(numpy.abs(nfft.fft(amp_cn)))**2/signal.size
    # print("gamma_max",gamma_max)

    # ---------计算特征值sigma_dp：非弱信号段零中心瞬时相位非线性分量的方差---------
    judge_A_n=amp_n>t_a
    sigma_dp=numpy.sqrt(numpy.sum(phi_NL[judge_A_n]**2)/numpy.sum(judge_A_n)-(numpy.sum(phi_NL[judge_A_n])/numpy.sum(judge_A_n))**2)

    # ---------基于以上特征值，识别调制方式ModulationType---------
    # modulationType:   =0，为常规调幅AM
    # modulationType:   =1，为DSBSC
    # modulationType:   =2，为SSB。在本实验中，为下边带LSB
    # modulationType:   =3，为FM
    modulationType=0
    if (sigma_dp<t_sigma_dp):
        modulationType=0 # 常规调幅AM
    else:
        if (abs(P)>=t_P):
            modulationType=2 # 为SSB, 本实验中，其为下边带LSB
        else:
            if (gamma_max<t_gamma_max):
                modulationType=3 # 为FM
            else:
                modulationType=1 # 为DSBSC

    # ---------计算信号频谱宽度---------
    if (modulationType==0):  # 若为常规调幅AM
        # 对于AM调制方式，由于其在载波频率处有一个大的冲激，其占据信号功率的绝大部分，
        # 因此，对于载频点在FFT中的位置fcn，必须计算完全准确，否则将给功率计算带来巨大误差，
        # 因此，为避免f0的少许误差使得fcn的值发生偏差，在上面已识别出是AM调制方式的前提下，
        # 下面重新计算fcn的位置
        # print("kc_origin",kc)
        kc=numpy.where(numpy.abs(signal_fft)==numpy.max(numpy.abs(signal_fft)))[0][0]
        # print("kc_changed", kc)
    # ---------计算总的功率---------
    TotalPower=numpy.sum(numpy.abs(signal_fft[0:round(signal.size/2)])**2) # 总的功率
    if (kc)<round(signal.size/2):
        TotalPower-=numpy.abs(signal_fft[kc])**2
    # print("kc=",kc)

    # ---------计算90%功率之和psum---------
    temp_value=numpy.abs(signal_fft[kc]) # 载频处的值在kc位置
    psum=0
    if (modulationType==0):
        psum=0 # 若是AM调制方式，psum不应包括载频处的功率，故psum初值为0
    else:
        psum= temp_value ** 2 # 初值为载频处的功率
    psum_stop=0
    # for i=kc:-1:1 # 从载频处开始计算功率之和
    for _ in range(kc-1,0,-1):
        if (modulationType==2): #若是单边带(下边带)调制，频谱不对称
            psum+=numpy.abs(signal_fft[_])**2
        else: #若是其他调制，频谱对称
            psum += numpy.abs(signal_fft[_])**2 #载频左边的功率
            psum += numpy.abs(signal_fft[kc + (kc - _)]) ** 2 # 载频右边的功率
        if (psum/TotalPower>=0.9):  # 信号带宽按总功率的90%来确定
            psum_stop=_
            break

    # ---------计算频谱宽度SpectrumWidth---------
    if (modulationType == 2):
        spectrumWidth = (kc - psum_stop) * f_Res
    else:
        spectrumWidth = (kc - psum_stop) * 2 * f_Res
    spectrumWidth = spectrumWidth / 1000 # 转化为单位：KHz

    # ---------显示结果---------
    modulationName=["AM","DSBSC","SSB","FM"]
    print("ModulationType={0}\nSpectrumWidth={1} KHz\nf0={2} Hz\nP={3}\ngamma_max={4}\nsigma_dp={5}\n".\
          format(modulationName[modulationType],spectrumWidth,f_signal,P,gamma_max,sigma_dp))


